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MONTE CARLO APPROXIMATIONS OF THE NEUMANN PROBLEM 



SYLVAIN MAIRE AND ETffiNNE TANRE 



Abstract. We introduce Monte Carlo methods to compute the solution of elliptic equa- 
tions with pure Neumann boundary conditions. We first prove that the solution obtained by 
the stochastic representation has a zero mean value with respect to the invariant measure of 
the stochastic process associated to the equation. Pointwise approximations are computed 
?H by means of standard and new simulation schemes especially devised for local time ap- 

^2 proximation on the boundary of the domain. Global approximations are computed thanks 

^>- to a stochastic spectral formulation taking into account the property of zero mean value of 

the solution. This stochastic formulation is asymptotically perfect in terms of conditioning. 

Cn Numerical examples are given on the Laplace operator on a square domain with both pure 

Cn Neumann and mixed Dirichlet-Neumann boundary conditions. 

r-| 1. Introduction 

-(— > 
C^ The aim of this paper is to compute pointwise and global numerical approximations of 

^ the solutions of pure Neumann problems for elliptic equations by means of Monte Carlo 

'~~' procedures. The pointwise solutions will be obtained via Feynman-Kac representations 

I and the global ones by stochastic spectral formulations inspired from the methods devel- 

^ oped in our previous works ifTSl [T9l . 

^^ We first consider the Neumann problem for the Laplace operator in a domain D dM.'' 

^~~^ which writes 

"^^ I ^Au{x) = -/(jc) for x^{xi,...,Xd)eD 

o 

^N where ^ stands for the incoming normal derivative of u on the boundary of the domain D. 

. . Because the boundary conditions are of Neumann type everywhere on the boundary, this 

J> equation has a solution up to an additive constant whenever an additional compatibility 

'k> condition on / and g is satisfied (see lilOJ ). This compatibility condition writes 

Cd [ fix)dx+ f giy)dy = 0. 

Jd JdD 

In the case of the Poisson equation with Dirichlet boundary conditions, the probabilistic 
representations of the solutions are based on Feynman-Kac formulae which are stopped at 
the first hitting time of the boundary by a Brownian motion. In the case of mixed boundary 
conditions, that is of Neumann type on a part of the boundary and of Dirichlet type on 
the other part, the (reflected) Brownian motion is still stopped at the first hitting time of 
the Dirichlet boundary. In the case of pure Neumann boundary conditions the Brownian 
motion never stops and the probabilistic representation introduced by Brosamler ||6l in the 
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case / = is the limit in time of the Feynman-Kac representation of the solution of the 
relative Cauchy problem. We have 

cT 



u{x) = lim E, ( f g{B,)dVs 



where (B? , f > 0) stands for the reflected Brownian motion in D and {Vi,t > 0) is its associ- 
ated local time on the boundary. The proof of this representation is based on probabilistic 
potential theory combined with a representation of additive functionals. This representa- 
tion has been extended by Bencherif-Madani and Pardoux |2l to more general advection- 
diffusion equations 

Lu{x) = -/(x), 

with boundary conditions -J^ = g{y), where 

, 1 ^ , .d^u{x) d, . .duix) 

z,/=l ^ J 1=1 ' 

and -j^ — ja'V{.).n stands for the conormal derivative. Some regularity assumptions are 
required on the domain D, the functions / and g, the coefficients a,- ^ and bj and also 
uniform eUipticity conditions and boundness for the symmetric matrix a. We have 

(1.1) u{x) = lim u{x,T) ^ lim E;^ ( / f{X,)ds+ / g{X,)dvA 

where {Xf.t > 0) is the reflected diffusion process and u{x,T) the solution of the cauchy 
problem associated to the operator L. The compatibility condition becomes 



f{x)p{x)dx+ g{y)piy)dy^O 

D JdD 

where p(x) is the density of the invariant measure associated to {Xi,t > 0). The proof of 
this more general representation is now based on exponential ergodicity. 

We intend to use the previous stochastic representation for the numerical computation 
of the Neumann problem by means of Monte Carlo procedures. We need to overcome 
several problems. To compute the poinwise solution, we have to choose a time Tq to stop 
the trajectories, deal efficiently with the approximation of the Neumann boundary condi- 
tions (which involves local time approximation) and understand which additive constant is 
obtained using this representation. To obtain the global solution, we need to choose and to 
build carefully the basis functions involved in the stochastic spectral formulations in order 
to obtain a unique and well characterized solution. The rest of the paper is organized as 
follows. 

In section l2] we first prove that the representation of Bencherif-Madani and Pardoux 
||2j naturally leads to a solution u{x) such that its mean value with respect to the invariant 
measure is zero. For a practical computation, we need to replace the limit in representation 
( |l.l| l by a finite time Tq. This introduces a bias that will be link to the second eigenvalue of 
the operator L with pure Neumann boundary conditions. The question of the variance of 
our estimators as a function of T is also studied. We remark that the variance has a linear 
growth as a function of T. 

In section [3] we describe a general algorithm to approximate the solution of general 
elliptic Neumann problems by means of the Euler scheme coupled with a local time ap- 
proximation method developed in lfT6l . 
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Section |4] is devoted to walk on spheres algorithms which can be used in the special 
case of the Poisson equation. We describe how to adapt them to a Cauchy problem with a 
finite horizon Tq. Furthermore, two new schemes with an increased order of convergence 
are introduced to deal with inhomogeneous Neumann boundary conditions. 

In section [5] we give some numerical results on mixed Dirichlet-Neumann problems 
and on pure Neumann problems with different degrees of difficulty. For the pure Neu- 
mann problem, the scheme that is used to compute the representation modifies the additive 
constant from which the solution depends. 

Finally in section l6] we develop stochastic spectral methods in order to obtain a very 
accurate and global approximation of the solution. In addition to the general methodology 
developed in ifTSl [T9l . we need to develop centering procedures for our approximation 
basis in order to project our solution in the space with zero mean value with respect to 
p{x). Numerical examples show that our method is very efficient in terms of accuracy and 
conditioning in the both cases where the invariant measure is known analytically or only 
approximated. 

2. Some properties of the stochastic representation 

For theoretical purpose, representation ( |l.l[ l is useful. For a numerical purpose, we need 
to understand which solution is effectively computed (not only up to an additive constant) 
and the consequences of replacing the limit by a fixed time Tq. We shall prove in the 
following that the solution u{x, T) has asymptotically a mean value zero with respect to the 
invariant measure and that its variance is increasing mainly linearly as a function of T. 

We consider a general diffusion operator L verifying the properties mentioned in the 
introduction |,2J. We denote by {Xt,t > 0) its associated reflected diffusion process 

Xr^x+ f b{X,)ds+ I a{X,)dW,-]- f riXs)dV, 
Jo Jo ^ Jo 

where the matrix a is such that GG^ — a, {Wi,t >Q) isad dimensional Brownian motion, 
{Vt,t > 0) is the local time on the boundary and y is in the conormal direction 

7(x) = an[x). 

We then define the solution u{x,T) of the Cauchy problem 

u{x,T)=Ejf f{Xs)ds+ f g{X,)dV, 
where / and g are bounded and verify the compatibility condition 



f{x)p{x)dx+ / g{y)p{y)dy^O. 

D JdD 

Following the arguments developed in |2|, there exists a bounded function /i such that 

u{x,T)=Ejf ,h{X,)ds 
with the new compatibility condition 

f\{x)p{x)dx — 
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2. 1 . Characterization of the solution. We can first notice that 



E.,(/i {X,))p{x)dx = E^ (/i (X,)) = / /i {x)p{x)dx = 

D Jd 

for any time s because jj. is the invariant law ofXs. Using Fubini theorem, this leads to 

/ m(x, T)p{x)dx = 
Jd 

for any time T. In the case of the Laplace operator, it has been proven in H] that the 
convergence in time of m(x, T) towards u{x) is uniform and thus the desired property holds. 
This uniform convergence relies on special properties of the Brownian motion that cannot 
be extended so easily to a general diffusion process. These properties are a bound on 
the transition density function p{t,x,y) of the reflecting Brownian motion and its spectral 
expansion |i7j. This spectral expansion writes 

Pit,x,y) = 7^ + E exp(-A„f)0„(x)0„(3;) 

I I n=l 

where rL is the uniform measure in D, (A„,n > 0) and (0„,n > 0) the eigenelements 
of —J A. The first term t^ corresponds to the first eigenvalue Xq — and ^o{x) — -4=. 
The remaining eigenvalues {X„,n > 1) are strictly positive. Theorem 2.4 of fT\ says that 
p{t,x,y) converges exponentially fast and uniformly to Jjj at a speed linked to the second 

eigenvalue — Ai of j A . This can guide us for the choice of the time Tq at which we should 
stop our trajectories. 

In the general case, we have for any T >Q, 



L 



u{x)p{x)dx 



u{x)p{x)dx — / u{x,T)p{x)dx 



< / {u{x) ~ u(x,T)) p{x)dx 

D 



thanks to the Cauchy-Schwartz inequality. In lemma 3 of ||2|, it is proven an exponential 
ergodicity for the solution w{x^t) = Ev(/i(X,)) of the parabolic problem with an initial 
condition f\ satisfying the centering condition, no source term and homogeneous Neumann 
boundary conditions. This is summarized in the inequality 

/ w(x,f)^/?(x)iix < exp(— rt) / fi{x)^p{x)dx 
Jd Jd 

where r is a positive constant. We can write 

{u{x) — u{x,T)) p{x)dx < / p{x)dx exp{-s)w{x,s)ds x / exp( — s)ds 
D Jd Jt ^ Jt ^ 

thanks again to Cauchy-Schwartz inequality and finally 

/ u{x)p{x)dx < -exp(— rr) / f\{x) p{x)dx 
Jd r Jo 

thanks to Fubini theorem and to the exponential ergodicity bound. This proves that J^u{x)p{x)dx 
0. 
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2.2. Estimation of the variance. In this section, we estimate the variance of /q fi{Bs)ds 
as a function of T in the simpHfied case of the Poisson equation in dimension one in the 
interval D = [a^b] with boundary conditions u (a) = u (b) = 0. We have 

2 / „T X 2 



Var[ / fi{B,)ds] = E, 







MSM" 



E, I fi{B,)ds 



where {Bs,s > 0) is a reflected Brownian motion in D. As 

the second term in the variance is obviously bounded because m is a continuous function in 
a bounded domain. For the first term, thanks to Ito-Tanaka formula we have 

f^ - f^ I - f^ - f^ I - 

u{Bt) = u{x)+ Lu{B,)ds+ u {B,)dB, ^ u{x) - fi{B,)ds+ u{B,)dB,. 
Jo Jo Jo Jo 

Thanks to the boundary conditions, the usual term involving the local time process is equal 
to zero. The above formula leads to 

MB,)ds] = {u{x)-u{BT)f+n u{B,)dB, 

+2 (u{x) -u(Bt)){ I u (Bs)dBs 
Using that the invariant measure is uniform in D, we have 

i2 1 







\\m¥.^[{u{x) - u{Bt))] 



\D\ 



(m(x) — m(>')) dy 



D 



which proves that this term is bounded and its limit identified. Then from the isometry 
property 



E, 







T ^ X 2- 

u'{B,)dB, 



E. 







u {Bs) ) ds 



and from the ergodic theorem 



E, 



M (B,)c/B., 



= rE, 



^ j (m'(B,)) ds 



^^iii /„("'«)'"■'■ 



Thanks now to the Cauchy-Schwartz inequality and to the previous bounds 



E;t {u{x) - u{Bt)) / M {B,)dB, 
which proves finally that asymptotically 



1 



2 \l/2 

<CVr(p^^(«(y)) dy 



C1-C2VT + C3T <Var{ / f{B,)ds)<Ci+C2VT+C3T 
Jo 

for some constants Ci , C2 and where C3 = t^t //) ( « (j) ) dy. The lower and upper bounds 

on the variance show that it increases essentially as a linear function of T (except if u{x) is 
constant). This proves that the time chosen to stop the trajectories is both crucial in terms 
of bias and in terms of variance. 
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3. EULER SCHEME APPROXIMATIONS 

We now present a scheme used to approximate the solution u inspired by the represen- 
tation ( |l.l[ l but where the limit is replaced by a fixed finite time Tq. This method relies 
on weak approximations. First of all, in a general setting, we have to approximate the 
reflected diffusion process {Xi,Q < t < Tq) associated to the infinitesimal generator L in 
the domain D C M''. The approximation of diffusion processes is usually done using the 
Euler scheme. The errors and the error expansions on the weak approximations using this 
scheme in the whole space have been first studied in |24j . The weak approximation of dif- 
fusions with absorbing (Dirichlet) and reflected (Neumann) boundary conditions have been 
treated respectively in |TT] and (5]. Given a time step 8, the approximation of {X,,t> 0) by 
a standard reflected Euler scheme {Xi^g,Q < k < To/5) can be described by the following 
procedure. 

(1) for all /=!..£/, 

;■ 

(2) ifl(A.+i)6 e A we set 1(^+1)3 ==l(i.+i)5, 

(3) else, if X(^+i)5 ^ D, we have to choose another position inside the domain D. 
The final position ^(^-+1)5 is the symmetrized of the Euler scheme ^(i^+i)^ in the 
conormal direction a.n (see [5J). 

To compute our representation, we first need to approximate the integral Jq " f{Xs)ds. This 
can be done using the rectangle method by the discrete sum 

S I f{x,s)- 

The approximation of the last term Jq° g{Xs)dVs is less classical (see for instance [16|| ). 
We introduce a delocalization parameter t, and compute the approximation 

8 L g{7l{Xk5))K^{X,s~K{X,5)) 
k=Q 

where Kf stands for instance for the Gaussian kernel 



We have also used the orthogonal projection n on the boundary dD. The delocalization 
parameter must be chosen carefully. We will observe in section [5] that it should not be too 
small unless the variance becomes very large. In this section, we will also consider a prob- 
lem with mixed boundary conditions. The treatment of the Dirichlet boundary conditions 
will rely on half-space approximations 1 12|. An additional stopping test of the trajectories 
based on a Brownian bridge enables to obtain an order 1 weak error in 5 on the part of the 
boundary with Dirichlet boundary conditions. 

4. Walk on spheres approximations 

4. 1 . Introduction. The Euler scheme can be used for the simulation of a wide class of 
stochastic processes linked to second order elliptic operators. If we deal with the Laplace 
operator or with divergence form operators with constant or piecewise constant diffusion 
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coefficients ifTTl . more efficient simulations of Brownian paths are available like walk on 
spheres (WOS) 1231 or walk on rectangles algorithms ||9l- 

In the case of Dirichlet boundary conditions, the stochastic representation of the solution 
implies a Brownian path up to the first time it hits the boundary dD of the domain. To solve 
the Neumann problem, we should use the same scheme as soon as the particle is away from 
the boundary, that is it has not reached the e-absorption layer Then, we replace the particle 
inside the domain and run again the same scheme until hitting one more time the boundary 
and so on. Some efficient ways to replace the particle inside the domain are described in 
section |43] 

Let us now describe the WOS method to compute the solution of the Laplace equation at 
a point X G D. The walk starts at x and jumps from one sphere to another until it reaches an 
e-absorption layer (= {y € D,d{y, dD) < e}). The spheres are built so that the jumps are as 
large as possible by taking the radius of the next sphere as the distance to the boundary dD. 
The next point is the exit first hitting point of the sphere by a Brownian motion started at 
the center. So, thanks to the isotropy of the Brownian motion, it has to be chosen uniformly 
on this sphere. We stop the walk at the first time the selected point on the sphere belongs 
to the e-absorption layer. The contribution of each walk to the solution is the value of the 
boundary term at a projection (generally the orthogonal one) of the current position on the 
boundary. Finally, the approximate solution is the average of the contributions of walks. 
In the case of the Poisson equation, it is also possible to compute the contribution of the 
source term during this walk fT5l using for each sphere S of radius r„ and center x„ a Green 
function conditioned by the exit point z writing 

(4.1) E,„ (^£m)dt\B,,=z^ = 'f l^z,y)fiy)dy. 

At least in dimension 2, the conditional density is known exactly and it is possible to 
sample easily from it. Thus we can obtain the contribution of the source term to the whole 
trajectory. 

4.2. WOS for the pure Neumann problem. Regardless the problem of making a good 
choice of a final time Tq to stop our trajectories, our task is to deal efficiently with the 
inhomogeneous Neumann boundary conditions while simulating accurately the elapsed 
time and position from the start of the trajectories. The description of the tools used when 
the walk hits the boundary is given in the next subsection. We need to adapt the WOS 



method to the case of a finite horizon Tq. First, formula (4.1 1 gives not the exit time T5 but 

^2 

only the average time -f spent in the sphere S of radius r„. Second, we want to approximate 
Jq^ f{Bt)dt knowing that B^^ = Z- This integral can be transformed writing J^^ f{Bt)dt — 
TsE^ [/ (Buzs)] where U is uniform on [0, 1] and independent of (B,,0 <t <ts) and E^ 
stands for the expectation with respect to this uniform r.v. U. This transformation called 
the one random point method has been used in our papers lITSl [T9l and also in ll22l in 
the context of financial mathematics. Its interest is to give a lot cheaper evaluation of the 
integral for only a small increase of the variance. 

Thanks to the isotropy of the Brownian Motion, we only need to simulate the couple 
(T5, jBjjt:^ ) knowing that B-^^ = (1,0) in the unit circle. To do this, we simulate an ab- 
sorbed Brownian motion starting at the center of the unit circle using the Euler scheme 
with a very small time step d and the half space approximation 1 12|. The walk stops after 
a random number M of steps at an exit point z — exp(/0 ) . The exit time Ts is approximated 
by M5, a point {x,y) is picked uniformly at random among the M points of the discretized 
trajectory and is rotated using a rotation of angle —0. We can thus obtain the empirical 
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joint law of (tjj ,B£/t5 ) knowing that B^^ ~ (liO). Samples from this empirical law are 
precomputed and stored in a large list. To sample from the couple (Ts^Bf/Ts \ we pick 
uniformly at random a couple (exit time, position) in the list and an angle a uniformly in 
[0, 2n\ to rotate the position of the previous couple. 

Dealing with the pure Neumann problem, we have to take care to another difference 
with the resolution of the Dirichlet problem: at the last step, the final time To is before the 
exit time of the sphere. Assuming that time is reinitialized to zero before the last sphere, 
we have to compute integrals of the form Jq ' f{Bs)ds knowing that B^^ = z and Ti < Ts- 
So, we also need to store Q full trajectories (time, position) before the exit time Tj to 
approximate the above integrals. 

4.3. Replacement after hitting the boundary. It remains to deal with the inhomoge- 
neous Neumann boundary conditions. When the process reaches the e-absorption layer it 
is projected on the boundary of D at a point {x,y). The standard way to replace the pro- 
cess inside the domain is based on a finite differences approximations Ii20j . The idea is 
to use a normal approximation of the derivative at the boundary. In order to simplify the 
description, we assume for instance that the boundary is locally vertical that is it has the 
form {{{x,z),ymm < z< ymax)}- Let {x,y) be a point on this part of boundary, we have the 
order one approximation 

u{x,y)~u{x + h,y) 

^ ^8{^,y) 

which leads to u{x,y) ~ hg{x,y)+u{x + h,y). This simply means in terms of randomization 
that hg{x,y) is added to the score of the walk and that the process starts again in the domain 
at position {x + h,y). This approach does not give the elapsed time from position (x,y) to 
position {x + h,y) and neither takes into account that we deal with the Laplace operator 

4.3.1. Kinetic approximation. A new scheme based on a kinetic approximation at the 
boundary has been introduced in [17| for the simulation of multidimensional diffusions 
in a media where the diffusion coefficient present some discontinuities. The homogeneous 
Neumann boundary conditions were also treated using this new scheme on a hard test case 
taken from the couplex exercices. This scheme is based on a small parameter approxima- 
tion of the diffusion operator by a neutron transport operator Another application was the 
Monte Carlo solution of the Poisson-Boltzmann equation in molecular dynamics ID. All 
the numerical tests were very satisfactory compared to the finite differences approach. The 
new scheme is proven to have an increased order of convergence on the test case of a single 
sphere for the Poisson-Boltzmann equation. 

We shall now describe this kinetic approximation in dimension 2 and how to adapt it to 
the treatment of inhomogeneous Neumann boundary conditions. We assume for the sake 
of simplicity that the process has hitten from the right a vertical boundary with a Neumann 
condition at a point {x,y) which is infinitely away from the other boundaries. We introduce 
a small parameter h in order to approximate the Laplace operator by the transport operator 
We also pick a collision time tc according to an exponential law of parameter 1 and a 
velocity (v^ = cos{9),Vy — sin(0)) uniformly on the half unit circle (— | < 9 < |). The 
new position is {x + hvJi.,y + h'^'ytc) and the elapsed time is h^t^. To deal with Neumann 
boundary conditions, we can write the Taylor expansion 

u{x + hvj„y + hvyt,) = u{x,y) +h\lfi{x,y) +h^\(f2{x,y) +0{h^) 

with 

f \ / N du(x,y) 

Vi [x,y) = ~vjc8{x,y) + Vytc — ^ — 
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and 

Taking its mean value, we obtain 

u{x,y)^muix + hv^tc,y + hvyt,)) + -^^^^h + fix,y)h^ + 0{lr') 

n 

because 

E(vv)=0, E(v.,)=/' cos{e)de, E(f2)=2. 

Our scheme is the randomization of the previous formula: the motion continues at point 

{x + hvJc,y + hvytc) and the quantity 

'^^h + fix,y)h^ 
n 

is added to the current score of the walk and h^ to the total time. The error is a 0{h^) 
each time the boundary is hitten during the walk if we do not take into account the other 
boundaries. If {x^y) is close to another part of the boundary, the point 

{x + hvJc,y + hvytc) 

may be outside the domain. In this case, we reduce h iteratively by a factor 2 until it is 
inside. In the case of a general boundary, the choice of the velocity law must be obviously 
adapted to the form of the boundary to ensure that the process reenters in D. 

4.3.2. Order three finite differences approximation. Some more sophisticated finite differ- 
ences schemes can be also be considered. Given a step h, we approximate the Laplace 
operator and the normal derivative at the point {x,y) of the boundary using order two finite 
differences. For the Laplace operator, we have thanks to the so-called diamond scheme 

1 -u{x-h,y)-u{x,y-h)-u{x,y + h)-u{x + h,y)+4u{x,y) 
--AM(x,y)=/(x,y)~ ^^ ■ 

at the order two and for the centered normal derivative 

du^ ^ u{x + h,y)-u{x-h,y) ^ ^(^^^) ^o(/^2)_ 
an 2h 

By combining these two equations and getting rid of the fictitious value u{x — h,y), we 
obtain 

2u{x + h,y)+u{x,y + h) + u{x,y-h) h^ \ , , r \ , r,tu'i\ 

u{x,y)^ -^ h— /(x,y)+/ig(x,y) + 0(/!-'). 

This formula can be used for Monte Carlo simulations. Each time the Brownian motion 
hits the boundary at a point (ji:,y),the quantity 

h^ 
-i^f{^.y)+hg{x,y) 

is added to the current score and the motion continues at one of the positions 

{x + h,y),{x,y + h),{x,y-h) 

with the discrete probability law (5,5,5). Some algorithms of a similar type for a walk 
on a grid have been introduced in |21| and are called sliding on the boundary methods 
because some of the possible points of replacement are on the boundary. Obviously if the 
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boundary is not a straight line, points {x,y + h) and (x,}? — /i) may lie outside the domain. 
If it happens, one needs to project the point on the boundary which induces an additional 
error. 

We introduce now a new similar scheme for a general boundary which avoids this pro- 
jection. Instead of writing the previous equations at point {x,y), we write them at point 
{x + h,y). We obtain 

—u(x,y) — u(x + h,y — h) — u(x + h,y + h) — u(x + 2h,y)+4-u{x + h,y) 
2h^ -/(- + /^,3') 

for the Laplace operator. For the normal derivatives, we write the two Taylor expansions 

/ , N / X , du(x,y) h^ d^u(x,y) ^,,^, 
uix + h,y) = u{x,y)+h \'^' + — ^^ + 0{h^) 
ox 2 dx-^ 



and 



u{x + 2h,y) = K^,}') +2/2^^^ +2/^'^^^ +0(/23) 

ax ax-'- 



which leads to 



du{x,y) -3u{x,y)-u{x + 2h,y)+Au{x + h,y) . 

g{x,y) = -^ — = ^7 + 0{h ). 

ox 2n 

We finally obtain 

u{x,y) ^ g{x,y)h + f(x + h,y)h^ + ^ '^ >-^ '^ >- + 0{h^) 

which leads to the following Monte Carlo representation. The motion is replaced equiprob- 
ably at point {x~\-h,y + h) or at point {x + h,y — h), 

gix,y)h+f{x + h,y)h^ 

is added to the total score and h^ to the total time. This time, the motion is more likely to 
be in the domain. If not, we reduce h iteratively by a factor 2 until it is, as we did for the 
kinetic approximation. 

5. Numerical results 

We shall now test our schemes on two Poisson equations in the square D = [—1,1] with 
increasing levels of difficulty. In the first one, we consider a mixed Dirichlet-Neumann 
where there is no problem of uniqueness for the solution and where the trajectories stop 
naturally. This enables a first comparison between the different schemes and validate the 
new ones introduced in section]?] In the second one, we study the pure Neumann problem 
for which additional difficulties arise, that we will try to overcome, especially the problems 
of the uniqueness of the solution and the choice of the time Tq to stop the walks. 

5.1. Mixed Dirichlet-Neumann Poisson equation. Our first test case is the Poisson equa- 
tion in the square domain D = [—1,1]^ with boundary dD = dD\ U dD2 defined by 

:= 1(^,3'), b| = 1,-1 <x<l}U{(x,y),x = -l,-l<y<l}, 
:={(x,y),x=l,-l<y<l}. 

We have a Neumann boundary condition on dDy and a Dirichlet boundary condition on 

dD2. 

-\Au{x,y) =/(x,y) := -a^exp(a(x + 3')) for {x,y)<ED 
\f„{x,y) =gi(x,y):=±fexp(a(x + y)) for {x,y)edDi 

u{x,y) =g2{x,y):=e.x^{a{x + y)) for {x,y)edD2, 
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(a,M) 


exact 


e{du^i) 


^ 

v^ 


e(52,&) 




e(53,<^3) 


03 


(ai,Mi) 


1.306 


9.9E-4 


2.0E-3 


4.7E-3 


1.7E-2 


3.8E-4 


2.0E-3 


(a2,Mi) 


1.705 


1.8E-3 


4.6E-3 


4.6E-2 


4.2E-2 


2.1E-4 


4.7E-3 


(0:3, Ml) 


2.226 


l.lE-2 


8.9E-3 


1.5E-1 


8.9E-2 


2.4E-3 


9.0E-3 


(ai,M2) 


1 


1.6E-3 


3.7E-3 


6.1E-2 


3.6E-2 


3.2E-3 


3.7E-3 


(«2,M2) 


1 


9.1E-3 


8.2E-3 


1.2E-2 


8.4E-2 


9.0E-3 


8.3E-3 


(a3,M2) 


1 


2.6E-2 


1.5E-2 


6.4E-2 


1.7E-1 


1.8E-2 


1.5E-2 


(ai,M3) 


0.766 


l.lE-2 


4.1E-3 


2.6E-2 


4.3E-2 


l.OE-3 


4.2E-3 


(a2,M3) 


0.587 


1.6E-2 


8.9E-3 


4.7E-3 


9.6E-2 


4.7E-3 


9.0E-3 


(a3,M3) 


0.449 


1.3E-2 


1.6E-2 


7.3E-2 


1.9E-1 


l.lE-2 


1.7E-2 



Table 1. Errors in the approximation of the solution of the 
Dirichlet-Neumann Poisson problem in the domain D = [—1,1] 
an Euler scheme with parameters (5i,(^i) = (0.01,0.01) , (^2 
(0.001,0.000001) and (^3,^3) = (0.001,0.001). 



mixed 
^ with 

^^2) = 



where the sign of ^1 is negative for y = +1 and positive otherwise. The stochastic process 
associated to this equation is a standard Brownian motion {Bt,t > 0), reflected on dDi and 
killed on dD2- The solution u has the stochastic representation 



u{x,y) = E(^j) 



td 



fiBMs' 



td 



giiB,)dV,+g2{Bro) =exp{a{x + y)) 



where To is the first hitting time of dD2- The parameter a is introduced to obtain solu- 
tions with different variations and consequently different variances in the Feynman-Kac 
representations. We choose the values ai — ^ , a2 — ^ and a3 = 1 ranked by increasing 
degree of difficulty. We compute the solution using A^ = 50000 Monte Carlo simulations 
at some reference points Ml = (0.8,0), M2 = (0,0) andM3 == (—0.8,0) which are located 
at positions at different distances from dD2. We also compute the ratio -j^, where ff^ is 
the variance of the method which gives an estimation of the Monte Carlo error 



5.1.1. Euler Scheme approximations. We denote by e{8,£,) the absolute error on the so- 
lution using the reflected Euler scheme with the half space approximation on the Dirich- 
let side. This error is computed in table [T] with 3 different sets of parameters (5i,(^i) = 
(0.01,0.01) , (52,<^2) = (0.001,0.000001) and (5?, .^3) = (0.001,0.001). 

For the parameters (5i,i^i), we observe a good accuracy of at least 2.6E-2 on the so- 
lution at all the reference points with the various levels of difficulty. We note that the 
approximate solution is more accurate and the variance smaller at point Mi which is very 
close to the Dirichlet boundary. The approximations at the two other reference points are 
similar in terms of error and variance. The CPU times on a standard PC for the compu- 
tation of the solution for the 3 levels of difficulty simultaneously are about 2 seconds for 
Ml, 7 seconds for M2 and 10 seconds for M3. Moreover, variance and bias increase with 
the level of difficulty. To improve the accuracy of our method, we have first chosen to 
reduce drastically the regularisation parameter and by a factor 10 the discretization pa- 
rameter The corresponding results for parameters (^2,^2) were not satisfactory because 
the variance increases too much with £,2- This is especially true when a ^ a^ where the 
accuracy is no more than one digit. The results are a lot better with (^3,^3) with roughly 
the same variance than with (5i,^i) but with a smaller bias. For example, the accuracy is 



MONTE CARLO APPROXIMATIONS OF THE NEUMANN PROBLEM 



12 



(a,M) 


^3(/^l) 


Fi{h2) 


K{h) 


K{h2) 


Fiih) 


Fiiln) 


a 


(ai,Mi) 


9.2E-4 


2.4E-3 


4.6E-3 


1.6E-3 


l.lE-2 


4.4E-3 


2E-3 


(a2,Mi) 


6.2E-3 


5.3E-3 


2.0E-2 


1.3E-3 


3.6E-2 


1.9E-2 


4E-3 


(aa.Mi) 


1.4E-2 


7.7E-3 


5.5E-2 


3.5E-3 


9.2E-2 


4.9E-2 


8E-3 


(ai.Mz) 


l.lE-2 


2.3E-3 


3.7E-3 


2.8E-3 


5.5E-2 


2.4E-2 


4E-3 


(a2,M2) 


1.3E-2 


1.8E-3 


3.9E-2 


3.3E-3 


1.2E-1 


8.5E-2 


8E-3 


(«3,M2) 


1.8E-2 


9.0E-3 


1.4E-1 


2.5E-2 


4.2E-1 


2.0E-1 


lE-2 


(ai,M3) 


1.7E-2 


l.lE-2 


2.5E-2 


1.2E-2 


6.3E-2 


4.2E-2 


4E-3 


("2,^3) 


1.2E-2 


l.OE-2 


4.1E-2 


l.lE-2 


2.7E-1 


1.4E-1 


8E-3 


(as, Mb) 


1.3E-2 


6.2E-3 


1.5E-2 


3.8E-3 


5.8E-1 


3.1E-1 


2E-2 



Table 2. Errors in the approximation of the solution of the mixed 
Dirichlet-Neumann Poisson problem in the domain D = [—1,1]^ with 
WOS approximations. The errors F^{h) and F2{h) are based on the finite 
differences method with scores respectively g{x,y)h + f{x + h,y)h^ and 
g{x,y)h at the boundary dD2- The error K{h) is computed thanks to the 
kinetic approximation, hi — 0.2 and h2 =0.1. 



3 times better at point Mi for the 3 values of a. The CPU times increase by a factor 10 
which corresponds to the reduction of 8. 

5.1.2. Walk on spheres approximations. In table |2j we compare three different methods 
to compute the absolute errors on the exact solution all relying on the walk on spheres 
method with absorption parameter e = 10^^ but with different ways to handle the Neumann 
boundary conditions. The first two errors Ft, [h] and F2 (h) are based on the finite differences 
method with scores respectively g{x,y)h+ f{x + h,y)h^ and g{x,y)h at the boundary dD2. 
This enables to emphasizes the differences between our new approach with the additional 
term f{x + h,y)h^ and the standard one. The last error K{h) is computed thanks to the 
kinetic approximation. The simulations are performed with two values hi = 0.2 and /12 = 
0. 1 of the parameter h. The exit time of the unit circle and the associated uniform position 
before absorption are pre-computed and stored in files of size 10^. 

The CPU times for the computation of the solution at points (Mi,M2,M3) with dis- 
cretization step hi are (4, 14,21) for the finite differences method, and (5,20,30) for the 
kinetic approximation (they are twice bigger when using /12). In terms of accuracy, we 
observe that the errors F2{h) are clearly a lot bigger than the ones obtained with the two 
new methods. The errors Fo,{h2) and K{h2) are similar and are furthermore comparable to 
e( 53,(^3) for a CPU time twice smaller. The error Fj,{lii) is similar to e(5i,^i) for a CPU 
time twice bigger. We can conclude that the two new schemes are very efficient especially 
when one desires an accurate solution. Nevertheless and unlike with Dirichlet boundary 
conditions, there is not such a big difference in terms of efficiency between Euler schemes 
and WOS methods. 

5.2. Pure Neumann Poisson Equation. 

5.2.1. Preliminary example. In this part, we illustrate our theoretical results of section|2] 
on the solution of a very simple Cauchy problem and of its related variance as a function 
of T. We consider the Poisson equation 



1 



Am = -f{x,y) = -2(3x2- 1)(/- 1)2-2(3/- l)(x2- 1) 
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Figure 5.1. Estimation of the solution of the pure Neumann problem 
at point {x,y) = (—0.5,-0.5) in function of the final time T 



for {x,y) & D = [—1 , l]^ and homogeneous Neumann boundary conditions ^ = on dD. 
The stochastic process associated to this equation is a standard reflected Brownian motion 
and its invariant measure is the uniform law in D. The exact solution with mean value zero 
is 



«(x,3')-(x2-l)2(/^l)2- 



64 

225' 



In this simple domain, we know that the second leading eigenvalue Ai of the operator 

. 2 

— jA with pure Neumann boundary condition is ^. This means that the convergence 
of the solution of the Cauchy problem towards the solution of the above equation is a 

2 

0(exp(— ^r)). We have also proven that the main part of the variance of our scheme 
increases linearly as CiT, where C3 = ^ /o |Vm| . Here, we have 



|V«|2 = (4x(x2 - 1)(/ _ 1)2)2 + (43;(/_ l)(x2 _ i)2)2 

and finally, C3 = 55573 — 0.99. We compute an approximate solution u{x,y,T) at point 
{x,y) — (-0.5,-0.5) and its related variance for values of T € [0.1,20]. The numerical 
method used is the Euler scheme with a small stepsize 5 — 0.001 and a huge number 
N = 5 X 10^ of simulations. The exact solution with mean value zero is 

m(-0.5,-0.5) ^^-^^ 0,03196. 



225 



In figure 5.1 we observe that m(0. 5,0.5,7) converges quickly to a constant (modulo some 



statistical variations) which is close to m(— 0.5, —0.5). 

Figure 15^ concerns the variance: we observe that it increases linearly as a function of 
T. If we compute the slope of the variance using for instance a linear regression on the 
approximate values at times 7) = 8 + / for 1 < / < 8, we obtain about 0.99. This means 
that the non linear part which behaves as C^s/T should be very small on this particular 
example. 

On this simple example, we have been able to confirm the results obtained in section [2] 
in the more general setting of a domain in dimension 2. However, we should now consider 
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Figure 5.2. Estimation of the variance of the approximation of the so- 
lution of the pure Neumann problem at point {x,y) = (—0.5,-0.5) in 
function of the final time T 



20 



inhomogeneous Neumann boundary conditions where both bias and variance may increase 
due to local time approximation. 

5.2.2. Numerical parameters and error criteria. Our main test case is the Poisson equa- 
tion in the square domain D — [—1,1]^ defined by 

1 , 

-- AM(jt:,3') = ~a e\p{a{x + y)) 

in D and with Neumann boundary conditions ^ |^ = ± " e\p{a{x + y)) on dD, where the 
sign is positive on the bottom and left sides of the boundary and negative on the right and 
top sides. The solution of this equation with zero mean value with respect to the invariant 
measure is hence 



u{x,y) —e\p{a{x+y)) — - 



exp {a{x + y)) dxdy 



-iJ-i 



that is 



u{x,y) =exp(a(x + 37)) — 



(exp(a)-exp(-a))^ 
4a2 



This solution is the one that we are likely to obtain numerically if we have a perfect sim- 
ulation of the reflected Brownian motion and a good choice of the time Tq when we stop 
the trajectories. Even though we have noticed on our test cases on the mixed problem 
that our numerical schemes are quite efficient, they are obviously not perfect and conse- 
quently introduce supplementary errors. In fact, regardless to the usual Monte Carlo and 
discretization error, we do not compute the solution of the equation with zero mean value 
with respect to the invariant measure but another one. Nevertheless, we know that our 
approximation Ua{x,y) should be of the form u{x,y) +a where a is a constant. In order 
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to check the quahty of our approximation and estimate this constant, we will compute the 
minimum of the weighted cost function 

Using a Gauss quadrature formula on the Tchebychef grid, J{a) can be approximated by 

1 

P2 



1 ^ ^ 

■^1 ('^) "= "^ L L Mxi,Xj) - u{xi,Xj) - a)^ 
i= 1.7=1 



using the Tchebychef points jc, = cos( ^ 'ip )■ The minimum of Ji is 

I P P 

" = "S2 E E MxnXj) - u{xi,Xj)) 

which indicates the bias with respect to the perfectly simulated solution while the value 
p = ^yJl{^) quantifies the adequacy to the model. The choice of the weighted cost func- 
tion is motivated by the use of bidimensional Tchebychef interpolation polynomials in the 
stochastic spectral methods of the next section. Indeed, they rely on the pointwise ap- 
proximations at the points of the Tchebychef grid. In practice, we choose P = 3 which 
is sufficient to have a good accuracy on the weighted integrals. Once a and p have 
been computed, we compute the approximate solution F{M) — Ua{x,y) and the error on 
the model B{M) — \ua{x,y) — u{x,y) — a\ for two sets of points M. The first one contains 
points M4 = (0,0), Ms = (—0.2,0.2), Me — (—0.8,0.8) and the second one contains points 
M^ = (0, 0.8), Mg = (0.2, 0.6), M9 = (0.4, 0.4) . Morevover all computations are performed 
using A^ ~ 50000 Monte Carlo simulations. 

5.2.3. Euler scheme approximations. First of all, we choose 7b = 10, which corresponds 
to a bias equal to exp(— 1.257r^) ~ 5 x 10^^ in all the following numerical tests. This 
value is small enough so that this bias is negligible with respect to the other errors and the 
variance not too large. 

The values of the quantities -j= depend weakly on the starting point and on the dis- 
cretization. They are approximatively equal to respectively 0.006, 0.013 and 0.03 for re- 
spectively ai = J, a2 = 5 and a^ — 1 which is about 1.5 times bigger than in the mixed 
Dirichlet-Neumann case. When (5i,(^i) = (0.01,0.01), a and p are equal respectively to 
-0.0059, -0.052, -0.130 and 0.0039, 0.027, 0.0498. When (82,^2) = (0.001,0.001), a 
and p are equal respectively to -0.0015,-0.0042, -0.016 and 0.0053,0.021, 0.042. 

We observe that the values of a and p are small which indicates a good adequacy to the 
model. Moreover, p and especially a are significantly closer to zero when the discretisation 
parameters decrease from (5i,^i) to (52,^2)- This shows that the approximate solution gets 
closer to u{x,y) as it has been proven in sectionl2] 

On this first set of points, we can see that for the same value of (a, 5,^) the direct 
estimations F{M4), F{Ms) and F{M(,) are close to each others. We observe also that the 
approximation model plays an important role. Indeed, the direct approximations are signif- 
icantly different from each other for the two sets of parameters but nevertheless the errors 
B{M) are quite small for both sets. The maximum absolute errors are 0.01, 0.038 and 0.031 
for respectively a\ , a2 and a^. 

The same conclusions hold for the second set of points. The maximum absolute errors 
are slightly bigger 0.02, 0.054 and 0.068 for respectively ai, a2 and a^. These maximum 
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F{M^) 


F{Ms) 


F{Me) 


B{Ma) 


B{Ms) 


B{Me) 


(«i,(5i,^i)) 


-0.045 


-0.042 


-0.042 


0.002 


0.001 


0.001 


(a2,(5i,^i)) 


-0.204 


-0.172 


-0.210 


0.005 


0.038 


0.001 


(a3,(5i,^i)) 


-0.480 


-0.506 


-0.508 


0.031 


0.005 


0.003 


(ai,(52,&)) 


-0.030 


-0.041 


-0.036 


0.01 


0.002 


0.003 


(«2,(52,fe)) 


-0.155 


-0.182 


-0.146 


0.006 


0.02 


0.016 


(«3,(52,&)) 


-0.374 


-0.365 


-0.395 


0.023 


0.032 


0.002 



Table 3 . Approximation F of the unbiased solution Ua of the pure Neu- 
mann problem with an Euler scheme with steps (5i,^i) = (0.01,0.01) 
and (§2,fe) = (0.001,0.001) at points M4 = (0,0), M5 = (-0.2,0.2), 
Mg = (—0.8,0.8) and the coiTesponding eiTors B{M) — 

\Ua{M) —u(M) — fl|. 





F{Mj) 


F{Mi) 


F{M<)) 


B(M7) 


B(M8) 


B{Mg) 


(ai,(5i,^i)) 


0.266 


0.255 


0.249 


0.004 


0.007 


0.013 


(a2,(5i,^i)) 


0.550 


0.524 


0.526 


0.054 


0.028 


0.031 


(a3,(5i,^i)) 


0.681 


0.715 


0.692 


0.033 


0.001 


0.022 


(ai,(52,^2)) 


0.249 


0.272 


0.271 


0.02 


0.005 


0.004 


(a2,(52,fe)) 


0.549 


0.537 


0.518 


0.006 


0.006 


0.025 


(a3,(52,fe)) 


0.832 


0.897 


0.790 


0.004 


0.068 


0.038 



Table 4. Approximation F of the unbiased solution m^ of the pure Neu- 
mann problem with an Euler scheme with steps (5i,^i) = (0.01,0.01) 
and (52, fe) = (0.001,0.001) at points M^ = (0,0.8), Mg = (0.2,0.6), 
Mg = (0.4,0.4) and the corresponding errors B{M) = 
\ua{M) — u(M) — fl|. 



errors are at least 2 or 3 times bigger than in the mixed Dirichlet-Neumann case. Further- 
more the CPU times are about twice larger in the pure Neumann Case. We can conclude 
that the pure Neumann problem is a lot harder to solve but that our algorithm still provides 
an acceptable accuracy for pointwise approximations. 



5.2.4. Walk on spheres approximations. We have noticed in section [J.1.2| that the two new 
methods to handle the boundary conditions have the same accuracy. We have chosen to 
use the order three finite differences in the following. Furthermore to make our simula- 
tions, we have precomputed and stored 100 positions of 100000 discretized trajectories as 
described in section |4~2| The time to open this file is negligible compared to the rest of the 
simulation times. They are about twice larger than in the mixed Dirichlet-Neumann case. 
The quantities -j= are still approximatively equal to respectively 0.006, 0.013 and 0.03 for 

respectively ai = j, a2 = | and a^ = I. For h\ = 0.1, a and p are equal respectively to 
-0.021, -0.106, -0.302 and 0.018, 0.053, 0.134. For hi = 0.05, a and p are equal re- 
spectively to -0.012, -0.060, -0.176 and 0.012, 0.034, 0.085. We observe that the value 
of p and especially a are larger than with the Euler Scheme method. 

Once again, the direct estimations F{Ma)^ F{Ms) and F{M^) are close to each others 
for a given value of the parameters {a,h). The maximum absolute errors are 0.006, 0.016 
and 0.049 for respectively ai, a2 and a^. This shows that even if the solution computed 
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F{Ma) 


F{Ms) 


F{M(,) 


B{Ma) 


B(M5) 


fi(M6) 


(ai,/!i) 


-0.055 


-0.059 


-0.062 


0.002 


0.001 


0.004 


(a2,/!i) 


-0.247 


-0.259 


-0.263 


0.016 


0.004 


0.001 


{a3,hi) 


-0.634 


-0.669 


-0.668 


0.049 


0.014 


0.015 


{oixM) 


-0.043 


-0.050 


-0.048 


0.006 


0.001 


0.001 


{a2,h2) 


-0.198 


-0.215 


-0.217 


0.019 


0.002 


0.001 


[a^M) 


-0.517 


-0.547 


-0.563 


0.04 


0.01 


0.005 



Table 5. Approximation F of the unbiased solution u^ of the pure 
Neumann problem with the WOS method at points Ma, = (0,0), M5 = 
(-0.2,0.2), Me = (-0.8,0.8) and the corresponding errors B{M) = 
\ua{M) — u{M) — fl|. The parameters are hi =0.1 and h2 = 0.05 



is further away from u{x,y), the accuracy is similar than the one obtained with the Euler 
scheme method. 





F{Mj) 


F{Ms) 


F(M9) 


B{Mj) 


BiMs) 


B(M9) 


(ai,/zi) 


0.231 


0.236 


0.240 


0.016 


0.011 


0.007 


(a2,/!i) 


0.399 


0.412 


0.424 


0.043 


0.03 


0.018 


(a3,/ii) 


0.455 


0.479 


0.503 


0.087 


0.063 


0.039 


(ai>2) 


0.246 


0.253 


0.253 


0.011 


0.003 


0.003 


(a2,/i2) 


0.466 


0.475 


0.474 


0.021 


0.012 


0.013 


(as, hi) 


0.627 


0.623 


0.626 


0.041 


0.045 


0.042 



Table 6. Approximation F of the unbiased solution m^ of the pure 
Neumann problem with the WOS method at points M^ — (0,0.8), 
Mg = (0.2,0.6), M9 = (0.4,0.4) and the corresponding errors B{M) = 
\u-a{M) — u{M) — a|. The parameters are h\ =0.1 and hi = 0.05 



The maximum absolute errors are bigger for this set of points especially for (a3,/ii). 
Nevertheless, the maximum absolute errors are 0.011, 0.021 and 0.045 for respectively 
ai, a2 and a^. We can conclude that we achieve a good accuracy on this pure Neu- 
mann problem but with an increased computational cost compared to the mixed Dirichlet- 
Neumann problem. 

6. Stochastic spectral methods 

6.1. Spectral formulation. In this section we describe how to adapt the stochastic spec- 
tral formulations introduced in 1 18 1 and studied in detail in 1 19| to the case of pure Neu- 
mann boundary conditions. These formulations are similar to usual spectral methods based 
on polynomial approximations |8 1 but they are built using relevant information, not neces- 
sarily at the collocation points, given by the Feynman-Kac formula. They are an extension 
of the sequential Monte Carlo algorithms for solving linear partial differential equations 
developed in lfT3l[T4l . These stochastic spectral formulations are asymptotically perfectly 
conditioned and quite easy to build. The case of mixed boundary conditions is nor de- 
scribed nor studied here because it is a straightforward extension of our previous works. 

When solving the pure Neumann problem using usual deterministic methods like finite 
elements, one also has to take into account very accurately the compatibility conditions and 
the non-uniqueness of the solution. Two approaches are usually used. The first one consists 
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in fixing the value of the solution at a specified node in order to avoid the resolution of a 
singular linear system. The second one leads to a singular system but it is solved using an 
iterative method like the conjuguate gradient for positive semi-definite linear system. In 
this second case, it is extremely important that the compatibility condition is verified at the 
discrete level which is obtained via discrete projectors. All these questions as well as an 
accurate study of the condition number are treated in |3|. 

Our stochastic formulation consists in computing a global linear approximation PpjU of 
the solution u using its values at some points x,- G D. This global approximation writes 

Af 
i=\ 

for some functions *I',(x) that are at least twice continuously differentiable and we assume 
that they verify the centering condition jj^^'i{x)p{x)dx — 0. Note that this last condition 
implies that P^u belongs to the space of functions that have a zero mean value with respect 
to p{x) and ensures the uniqueness of the solution. We also assume that for every point x;, 
we can approximate m(x,) via for instance a numerical approximation of the Feynman-Kac 
formula by 

SIx,) = ta,^g(zff") + f ft.;/(4f"). 
In practice this approximation is also such that 

lim m(x,) = m(x,) — / u{x)p{x)dx 
p,</,r^oo,5^o Jo 

which indicates that the solution that we compute numerically is close to the one with 
zero mean value with respect to p{x). The coefficients a,,y and j3,./ are positive weights. 
The p points z, ' ' are located on the boundary dD, the q points z^' ■ '' mD, 8 stands for the 
discretization parameter of the simulated reflected diffusion and T is the deterministic time 
when we decide to stop our random walk. We now let rN{x) — u(x) ~Pnu{x) and write the 
partial differential equation solved by r^ (x) . We have 

Lr^ = Lu — LP^u = / — LPnu 

in D with boundary conditions 

drN dPf^u 



and hence the approximation 




dria dria 



dPm(J:f)\ , f^ / ,,,5,,, 



which leads to the linear system Cu — d with 



S.T^ 



for i y^ k, 



=1 ^"^ ,=i 



;=i '^"« i=i 
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and 

7=1 7=1 

As we have done in ifTSl [T9ll . we can look at the asymptotic system we obtain when 
p,q,T ^oo and when 5^0. The term 

;=i "^"^ )=i 

is our Monte Carlo approximation at point jc, of the solution of the equation 

Lv = L'i'k 

with boundary conditions 

dv (9^,. 



dria dria 
on (9D that is ^k{xi) because j^^^(x)/?(x)t/jc = 0. As in our previous papers, this shows 
immediately that the matrix of the asymptotic system converges toward the identity matrix 
of size A^. This also means that the condition number of the system is naturally close to one 
even without additional preconditioning techniques like for instance the ones developed in 
123. 



6.2. Centering procedures. 

6.2.1. Exact centering. We now describe how this method works in practice. The main 
problem is that in general, usual linear approximations do not verify the centering condi- 
tions. We consider for the moment that we start with A^ + 1 Lagrange interpolation polyno- 
mials (pi at A^ + 1 points x,- G D. Such functions verify (pi{xk) — 5ik- The usual polynomial 
interpolation Qm+iu of degree A^ of a function u writes Qn+iu{x) — L^^' u{xi)(pi{x). Con- 
sidering the constant function ii = l,we obviously have L^' <Piix) = 1 and hence 

N+l 



v + l r 

E / (Pi{x)pix)dx^ 1, 
,= 1 JD 



= 1 JD 

which proves that this usual interpolation cannot verify the centering conditions. Never- 
theless, we can chOose an index /o such that 



9io{x)p{x)dx 

D 



max 

l<KAr+l 



^i{x)p{x)dx 

D 



^0 



for a sake of stability in the following approximation. The centering condition for Qn+\_u 
writes 

N+\ 



Y^ ( / (Pi{x)p{x)dx \u{xi)=0 



7=1 

which leads to 

Letting now ^,(x) = (pilx) — ji^^'y y (pi„(x), the new basis functions verify the center- 
Jo 'P'o^'^-'^^'^' 
ing conditions Jjy^i{x)p{x)dx — and still ^/(x^.) — 5, ,t. This centering procedure can be 

easily extended to general linear approximations not necessarily of interpolation type. It 
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is essentially the same approach than the projection method described in O for the finite 
elements method. 

6.2.2. Approximate centering. This solves our problem whenever the integrals J^ (pi{x)p{x)dx 
can be computed exactly. This may happen when the domain is simple and when the den- 
sity p{x) is known. This include for example the case of the Poisson equation in a hy- 
percube using an interpolation on a Tchebychef grid. In many other cases, these integrals 
need to be computed numerically. The density p{x) can be approximated by the law of the 
position y of M particles moving according to the reflected diffusion starting at any given 
point in D at a time Ti large enough. We obtain 



f 1 ^ 

/ (p,{x)p{x)dxC:^-Y,(p,iYj) 

Jd m . , 



and new basis functions which are defined by 



The coefficients of the asymptotic spectral matrix are 

j Q,; = 1 - Jj)^i{x)p{x)dx 

[cij =- jjj^i{x)p{x)dx for i^j. 

An easy computation shows that its eigenvalues are all but one equal to one. The remaining 
eigenvalue is 



N j- 

^N = l-Y. 'Vi{x)p{x)dx 

i=\JD 



and we have 

and the inequality 

|1-A^ 



E / (Pi{x)p{x)dx- f l^ ' ' Jx / (Pio{x)p{x)dx 
i=i\JD mLj=i%AYj)Jd j 



J M 



M 



i=i 



< 



/•AT J M AT 



1 



M 



(pi^{x)p{x)dx)-^Y.^'oiYj) 



(Pig{x)p{x)dx 



J^(Pi{x)p{x)dx 

Di=\ 



which proves that Xn converges to 1 when M -^ oo and that the condition number is once 
again asymptotically one. Note that we compute with this formulation an approximation 



Pnu{x) = ^m(x;)^;(x) 

of the solution with discrete integral equal to zero with respect to the particle approximation 
ofp{x). 
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6.3. Application to the Poisson equation. In this section, we describe the application 
of the stochastic spectral formulation on our main example of the Poisson equation with 
pure Neumann boundary conditions studied in section l5] The basis functions will rely on 
Tchebychef interpolation polynomials in dimension 2. The big advantage of this test case 
is that the centering procedure can be done exactly as we know that the invariant probability 
is the uniform law and because the integration domain is a square. This also enables us 
to compare the exact centering procedure and the approximate one where the numerical 
integration is done by means of a particle approximation of the invariant measure. 

6.3.1. Basis functions. Our spectral approximation is based on the standard tensorized 
interpolation of the solution on the Tchebychef grid which writes 

N N 
71=0 m=0 

where /„ is the Lagrange polynomial associated to z„ and m„ „, is the approximate value 
of the solution at the point (z„,z,„) where Zj — cos ( 2y(, , 2 ?r), < j < N. For the sake of 
simplicity, we choose A^ even so that point (0,0) belongs to the Tchebychef grid. Indeed, 
in this case the maximum of the integrals 

1 
l„{x)lm{y)dxdy 
-iJ-i 

is always attained for the Lagrange polynomials In{x) and lN{y) corresponding to this 

point. The function In {x)In (y) is removed from the basis functions and the (A^+ 1)^ — 1 

centered basis functions ^„ ,„(x,y) now write 

/_ 1 /_ 1 In {x)lm {y)dxdy 

'i'n.m{x,y) = ln{x)lm{y) - , , lN{x)lN{y) 

J_iJ_ilN[x)lN(y)dxdy ^ ^ 

for («,m) ^ {j,j). In the case of the approximate centering, the integrals above are re- 
placed by their particle approximations. 

6.3.2. Numerical results: exact centering. We present our results based on the Euler scheme 
approximation with two time discretization parameters 5i = 0.01,^2 = 0.001, a regular- 
ization parameter ^ = 0.001 and with two different number of simulations Mi = 1000 
and M2 — 5000. The trajectories are stopped at final time Tq = 10. For these four sets of 
parameters, we compute the maximum absolute error over the grid points 

erri{N) ~ max |m(z,-,z;) — m; J 

0<i,j<N ' 

and the condition number k{N) of the spectral matrix for Ni—2 and A^2 = 4 and a= \. 

We observe that the condition number K is decreasing as M increases and 5 decreases. 
The system is very well conditioned especially for the parameters (52, M2). As for usual 
standard spectral methods, the error is small and decreases with K and when A^ increases. 

6.3.3. Numerical results: approximate centering. The spectral method now requires the 
approximation of the invariant measure which should be done with Q simulations and a 
time step 5. The invariant measure is the uniform law in the square [—1,1]^. To study its 
impact on the spectral matrix, we have chosen to make its approximation using simply Q 



MONTE CARLO APPROXIMATIONS OF THE NEUMANN PROBLEM 



22 





err\{N\) 


K{N,) 


err\{N2) 


K{N2) 


(5i,Mi) 


6.1E-3 


11.4 


3.6E-5 


856 


(5i,M2) 


3.4E-3 


2.9 


1.2E-5 


121 


(52,Mi) 


4.2E-3 


5.1 


2.3E-5 


114 


(52,M2) 


l.lE-3 


2.3 


3.1E-6 


17 



Table 7. Global error err\ [N) — maxo<,j<Af |m(z;,Z;) — M;,;|of the ap- 
proximation of the solution of the pure Neumann problem with Euler 
scheme with time steps 5\ — 0.01,52 = 0.001, and a regularization pa- 
rameter ^ = 0.001 for for Ni—2 and A^2 = 4 and a = \. The value of 
condition number of the spectral matrix is k{N). 





erri{Ni) 


err2 {Ni ) 


K{Nl) 


err\{N2) 


err2{N2) 


K{N2) 


(5i,Mi,ei) 


1.9E-2 


2.4E-3 


10.2 


4.60E-2 


2.7E-5 


550 


{8uMi,Q2) 


5.8E-3 


5.1E-3 


26.3 


1.5E-3 


2.3E-4 


1187 


{Si,M2,Qi) 


1.7E-2 


1.6E-3 


5.2 


1.4E-2 


6.3E-5 


1130 


(5i,M2,e2) 


5.1E-3 


6.1E-4 


2.9 


3.1E-3 


l.OE-5 


644 


(52,Mi,ei) 


1.6E-2 


3.6E-3 


4.1 


1.5E-2 


5.4E-5 


309 


(52,Mi,e2) 


3.9E-3 


1.4E-3 


3.0 


7.4E-4 


1.4E-5 


252 


(52,M2,ei) 


6.3E-2 


9.4E-4 


2.8 


3.0E-2 


3.1E-6 


22 


(52,M2,e2) 


1.5E-3 


1.2E-3 


1.6 


3.3E-3 


4.2E-6 


16 



Table 8. Global approximation of the solution of the pure Neumann 
problem with an Euler scheme. The number of points are respectively 
N\=2 and N2 = 4. The time steps are 5i = 0.01,^2 = 0.001. The 
invariant measure is approximated with samples of sizes Qi = 100 and 
22 = 10000. 



samples {Xi,Yi) of this uniform law. We use two samples of different sizes Qi = 100 and 
Q2 — 10000 and introduce a new error criterion 



err2(N) — max 
0<ij<A' 



1 2 

u{zi,Zj) ~ Uij + 27; E «(^''^') 



as we have proven that we approximate the solution with discrete integral equal to zero 
with respect to the particle approximation of p{x). We nevertheless keep also the previous 
error criterion to study the impact of the particle approximation on the bias. 

The key observation is that we check that the solution effectively computed is the one 
with discrete integral equal to zero with respect to the particle approximation of p{x). Both 
accuracy and condition number have the same behaviour with respect to the parameters 5 
and M than with the exact centering. The condition number and erri decrease when the 
particle approximation is done with more points. 



7. Conclusion 

To compute Monte Carlo approximations of the solution of the Neumann problem for 
elliptic equations, we had to overcome several difficulties. First, we have characterized the 
solution of the Feynman-Kac representation introduced in ^ as the one with zero mean 
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value with respect to the invariant measure of its associated stochastic process. Then, we 
have proven that the variance increases mainly linearly as a function of time. 

We have introduced some new schemes to deal with the inhomogeneous Neumann 
boundary conditions. They were tested first on a pointwise approximations of mixed 
Dirichlet-Neumann problem where they show a good efficiency. The pointwise resolu- 
tion of the pure Neumann problem was a lot harder. Indeed, we had to choose the time Tq 
to stop our trajectories not too large because of the increase of the variance but also not 
too small because of the bias. We had also to understand that the solution computed nu- 
merically depends on the parameters of the numerical schemes and is not equal to the one 
with zero mean value with respect to the invariant measure. Taken all these difficulties into 
account, we have been able nevertheless to reach a reasonable accuracy on the approximate 
solutions. 

Concerning the global spectral approximation, we had to pay attention to the zero mean 
value property of the solution to chose our approximation basis. This has been achieved us- 
ing exact or approximate centering procedures very similar to the usual ones used in finite 
element methods. In both cases, the condition number of the spectral matrix was proven 
to be asymptotically one. The numerical experiments show that the stochastic spectral 
method is both very accurate and well-conditioned. 

The pointwise approximations of the pure Neumann problem are not completely satis- 
factory because the solutions obtained depend on the parameters of the numerical scheme. 
The choice of the time 7b to stop the trajectories is not straightforward for a general dif- 
fusion in a complex domain. For the WOS method, it requires furthermore to keep in 
memory discretisation of trajectories which is both costly and adds an error not so easy to 
quantify. It could be interesting to add a penalization term either in the source term or via 
robin boundary conditions to at least get rid of some of these drawbacks. 
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